---
title: "Politicized and Sensitive Issues"
output: html_document
---

Note: Chunks must be run in order.

# Setup

```{r}
rm(list=ls())

require(tidyverse)
require(broom)
require(stringr)
require(stargazer)
require(ggrepel)
require(pscl)

# load data
load("replication_data.RData")
```

## Reshape

Here we reshape the data to make it easier to model each issue independently.

```{r}
# save original issue names
issue_names <- dyad_level %>% 
  dplyr::select(`Asylum Seekers and Refugees`:`Women's Rights`) %>% 
  names()

# participating dyads only
by_issue_part <- dyad_level %>% 
  filter(participation == 1) %>%
  gather(issue, n_issue, `Asylum Seekers and Refugees`:`Women's Rights`) %>%
  dplyr::select(1:6, issue, n_issue, 7:length(.)) %>%
  group_by(issue) %>%
  nest() 
```

## Model Functions

Here we define some helpful functions for statistical modeling.

```{r}
# main regression
issue_model_part <- function(df) {
  mod <- lm(n_issue ~ ip_affinity_lagged + n_recs + mean_severity +
              hr_latentmean_lagged.x + hr_latentmean_lagged.y +
              same_region + upr_reviewed_session +
              factor(year) + factor(To_COW) + factor(From_COW), data = df)
  
  print("finished")
  
  # extract effect of geopolitical affinity
  geo_efx <- tidy(mod) %>% 
    filter(term == "ip_affinity_lagged") %>%
    dplyr::select(estimate, std.error)
  
  return(geo_efx)
}

# marginal effects function:
# extracts marginal effects and confidence intervals
# for geopolitical affinity
efx <- function(models){
  modelcoef = round(as.numeric(models$estimate), 4)
  modelse =round(as.numeric(models$std.error), 4)
  ylo = modelcoef - qnorm((1-0.95)/2)*(modelse)
  yhi = modelcoef + qnorm((1-0.95)/2)*(modelse)
  names = models$issue
  sign = ifelse(ylo < 0, "red", ifelse(yhi > 0, "blue", "gray"))
  marg.efx = data.frame(names, modelcoef, modelse, ylo, yhi, sign)
  
  return(marg.efx)
}  
```


## Plotting Functions

Here we define functions that we'll use for plotting.

```{r}
# plotting function - politicized issues
coefplot.politicized <- function(dfplot){

  #define plot
  p = ggplot(dfplot, aes(x = fct_reorder(names, modelcoef, .desc = TRUE), 
                         y = modelcoef, ymin = ylo, ymax = yhi)) + 
    geom_pointrange(colour = dfplot$sign) + 
    geom_hline(yintercept=0, color = "black") +
    theme_bw() +
    ylab("Marginal Effect of Geopolitical Affinity") +
    xlab("Issue") +
    scale_y_continuous(breaks = scales::pretty_breaks(n = 10)) +
    coord_flip()
  return(p)
}

# plotting function - sensitivity issues
coefplot.sensitivity <- function(models){
  modelcoef = as.numeric(models$estimate)
  modelse = as.numeric(models$std.error)
  ylo = modelcoef - qnorm((1-0.95)/2)*(modelse)
  yhi = modelcoef + qnorm((1-0.95)/2)*(modelse)
  names = models$issue
  dfplot = data.frame(names, modelcoef, modelse, ylo, yhi)
  
  #define plot
  p = ggplot(dfplot, aes(x = fct_reorder(names, modelcoef, .desc = TRUE), y=modelcoef, ymin=ylo, ymax=yhi)) + geom_pointrange(colour=ifelse(ylo < 0, "red", ifelse(yhi > 0, "blue", "gray"))) + 
    geom_hline(yintercept=0, color = "black") +
    theme_bw() +
    ylab("Marginal Effect of Issue on Likelihood of Support") +
    xlab("Issue") +
    scale_y_continuous(breaks = scales::pretty_breaks(n = 10)) +
    coord_flip()
  return(p)
}
```

# Figure 2: Effect of Geopolitical Affinity on Recommendation Counts, by Issue

```{r}
# model all issues
by_issue_mods <- by_issue_part %>% 
  mutate(model = map(data, issue_model_part)) %>%
  dplyr::select(-data) %>%
  unnest(model)

# get effects
effects <- efx(by_issue_mods)

# save for later
write.csv(effects, "Results/Figure_2/F2_effects.csv", row.names = F)

# plot
coefplot.politicized(effects)
ggsave("Results/Figure_2/F2_plot.pdf", width = 8, height = 8)
```

# Table A6: Effect of Geopolitical Affinity on Recommendation Counts, by Issue

Regression Tables for Figure 2

```{r}
# re-define modeling function to return full models
issue_model_untidy <- function(df) {
  print("finished")
  mod <- lm(n_issue ~ ip_affinity_lagged + n_recs + mean_severity +
              hr_latentmean_lagged.x + hr_latentmean_lagged.y +
              same_region + upr_reviewed_session +
              factor(year) + factor(To_COW) + factor(From_COW), data = df)
  return(mod)
}

# model all issues
by_issue_mods <- by_issue_part %>% 
  mutate(model = map(data, issue_model_untidy)) %>%
  dplyr::select(-data) 

# write tables
for(i in c(1:6)){
  print(str_c("writing: ", i))
  stargazer(by_issue_mods$model[(((i-1)*9)+1):(i*9)], 
            out = str_c("Results/Table_A6/Table_A6-", i, ".html"), 
            type = "html", 
            style = "apsr", 
            title = "", 
            # se = list(logit1.se, logit2.se, logit3.se), 
            notes="Reviewer, Target, and Year Fixed Effects Omitted from the Table. Standard errors in parentheses.",  
            dep.var.labels = "Number of Recommendations", 
            column.labels = by_issue_mods$issue[(((i-1)*9)+1):(i*9)],
            model.numbers = F,
            keep = 1:8,
            keep.stat = c("n"),
            covariate.labels=c("Geopolitical Affinity", 
                               "Total Recs", 
                               "Mean Severity", 
                               "Physical Integrity Rights (Reviewer)",
                               "Physical Integrity Rights (Target)",
                               "Same Region",
                               "UPR Review (Reviewer)"), 
            star.cutoffs = c(0.10, 0.05, .01))
}
```

# Figure 3: Effect of Geopolitical Affinity on Recommendation Counts, by Theme 

```{r}
# reshape using thematic cluster 
by_issue_cluster <- dyad_level %>% 
  filter(participation == 1) %>%
  gather(issue, n_issue, starts_with("cluster_")) %>%
  select(1:6, issue, n_issue, 7:length(.)) %>%
  group_by(issue) %>%
  nest() 

# model all issues
by_issue_mods <- by_issue_cluster %>% 
  mutate(model = map(data, issue_model_part)) %>%
  select(-data) %>%
  unnest(model)

by_issue_mods$issue <- c("Civil-Political Rights", "Governance and Public Services", "Migration", "Physical Integrity Rights", "Racial, Ethnic, and Religious Minorities", "Socio-Economic Rights", "Protection of Vulnerable Populations")

# get effects
effects <- efx(by_issue_mods)

# plot
coefplot.politicized(effects) + 
  theme(text = element_text(size=13))

# save
write.csv(by_issue_mods, "Results/Figure_3/F3_effects.csv", row.names = F)
ggsave("Results/Figure_3/F3_plot.pdf", width = 8, height = 4)
```

# Figure A1: Zero-Inflated Poisson: Count Coefficients for Geopolitical Affinity, by Issue

```{r}
# reshape using all possible dyads
by_issue_all <- dyad_level %>%
  gather(issue, n_issue, `Asylum Seekers and Refugees`:`Women's Rights`) %>%
  mutate(n_issue = replace_na(n_issue, 0), # replace NAs with 0
         n_recs = replace_na(n_recs, 0)) %>%
  dplyr::select(1:6, issue, n_issue, 7:length(.)) %>%
  group_by(issue) %>%
  nest()

# define zeroinfl model
issue_model_zeroinfl <- function(df) {
  mod <- zeroinfl(n_issue ~ ip_affinity_lagged + n_recs + mean_severity +
              hr_latentmean_lagged.x + hr_latentmean_lagged.y +
              same_region + upr_reviewed_session | 
              HRC_member.x + upr_reviewed_session + hr_latentmean_lagged.x + same_region,
              dist = "poisson",
              data = df)
  
  print("finished")
  
  # tidy
  x <- data.frame(type = c("count"),
                  estimate = c(coef(mod, "count")["ip_affinity_lagged"]),
                  ylo = c(confint(mod)["count_ip_affinity_lagged", 1]),
                  yhi = c(confint(mod)["count_ip_affinity_lagged", 2]))

  return(x)
}

# model all issues
by_issue_mods <- by_issue_all %>% 
  mutate(model = map(data, issue_model_zeroinfl)) %>%
  dplyr::select(-data) %>%
  unnest(model) 

# get counts
counts <- by_issue_mods %>%
  filter(type == "count") %>%
  ungroup() %>%
  mutate(names = issue_names,
         modelcoef = estimate,
         sign = ifelse(yhi < 0, "red", ifelse(ylo > 0, "blue", "gray"))) %>%
  dplyr::select(names, type, modelcoef, ylo, yhi, sign )

# plot
coefplot.politicized(counts) +
  ylab("Estimated Coefficient of Geopolitical Affinity")

# save
ggsave("Results/Figure_A1/FA1_plot.pdf", width = 8, height = 8)
write.csv(counts, "Results/Figure_A1/FA1_effects.csv", row.names = F)
```

# Figure A2: HRC members only

```{r}
# reshape using HRC members only
by_issue_hrc <- dyad_level %>% 
  filter(HRC_member.x == 1) %>%
  gather(issue, n_issue, `Asylum Seekers and Refugees`:`Women's Rights`) %>%
  dplyr::select(1:6, issue, n_issue, 7:length(.)) %>%
  group_by(issue) %>%
  nest() 

# model all issues
by_issue_mods <- by_issue_hrc %>% 
  mutate(model = map(data, issue_model_part)) %>%
  select(-data) %>%
  unnest(model)

# get effects
effects <- efx(by_issue_mods)

# plot
coefplot.politicized(effects)

# save
write.csv(effects, "Results/Figure_A2/FA2_effects.csv", row.names = F)
ggsave("Results/Figure_A2/FA2_plot.pdf", width = 8, height = 8)
```

# Figure A3: Effect of Geopolitical Affinity on Recommendation Counts, by Issue (Severity Level 3 Only)

```{r}
# reshape using severity 3 only
by_issue_3 <- dyad_level_S3 %>%
  filter(participation == 1) %>%
  gather(issue, n_issue, `Asylum Seekers and Refugees`:`Women's Rights`) %>%
  mutate(n_issue = replace_na(n_issue, 0),
         n_recs = replace_na(n_recs, 0)) %>%
  select(1:6, issue, n_issue, 7:length(.)) %>%
  group_by(issue) %>%
  nest()

# model all issues
by_issue_mods <- by_issue_3 %>% 
  mutate(model = map(data, issue_model_part)) %>%
  select(-data) %>%
  unnest(model)

# get effects
effects <- efx(by_issue_mods)

# plot
coefplot.politicized(effects)

# save
write.csv(effects, "Results/Figure_A3/FA3_effects.csv", row.names = F)
ggsave("Results/Figure_A3/FA3_plot.pdf", width = 8, height = 8)
```

# Figure A4: Effect of Geopolitical Affinity on Recommendation Counts, by Issue (Formal Alliance)


```{r}
# define alliance model
issue_model_alliance <- function(df) {
  mod <- lm(n_issue ~ alliance_lagged + n_recs + mean_severity +
              hr_latentmean_lagged.x + hr_latentmean_lagged.y +
              same_region + upr_reviewed_session +
              factor(year) + factor(To_COW) + factor(From_COW), data = df)
  
  print("finished")
  
  # extract effect of geopolitical affinity
  geo_efx <- tidy(mod) %>% 
    filter(term == "alliance_lagged") %>%
    dplyr::select(estimate, std.error)
  
  return(geo_efx)
}

# model all issues
by_issue_mods <- by_issue_part %>% 
  mutate(model = map(data, issue_model_alliance)) %>%
  select(-data) %>%
  unnest(model)

# get effects
effects <- efx(by_issue_mods)

# plot
coefplot.politicized(effects)

# save
write.csv(effects, "Results/Figure_A4/FA4_effects.csv", row.names = F)
ggsave("Results/Figure_A4/FA4_plot.pdf", width = 8, height = 8)
```


# Figure A5: Hurdle Count Models 

```{r}
# define hurdle model
issue_model_hurdle <- function(df) {
  mod <- hurdle(n_issue ~ ip_affinity_lagged + n_recs + mean_severity +
              hr_latentmean_lagged.x + hr_latentmean_lagged.y +
              same_region + upr_reviewed_session,
              dist = "poisson",
              zero.dist = "binomial",
              data = df)
  
  print("finished")
  
  # tidy
  x <- data.frame(type = c("count", "zero"),
                  estimate = c(coef(mod, "count")["ip_affinity_lagged"], 
                               coef(mod, "zero")["ip_affinity_lagged"]),
                  ylo = c(confint(mod)["count_ip_affinity_lagged", 1],
                          confint(mod)["zero_ip_affinity_lagged", 1]),
                  yhi = c(confint(mod)["count_ip_affinity_lagged", 2],
                          confint(mod)["zero_ip_affinity_lagged", 2]))

  return(x)
}

# model all issues
by_issue_mods <- by_issue_part %>% 
  mutate(model = map(data, issue_model_hurdle)) %>%
  dplyr::select(-data) %>%
  unnest(model) 

# get effects for zeros
zeros <- by_issue_mods %>%
  filter(type == "zero") %>%
  ungroup() %>%
  mutate(names = issue_names,
         modelcoef = estimate,
         sign = ifelse(yhi < 0, "red", ifelse(ylo > 0, "blue", "gray"))) %>%
  dplyr::select(names, type, modelcoef, ylo, yhi, sign )

# plot
coefplot.politicized(zeros) +
  ylab("Estimated Coefficient of Geopolitical Affinity")

# save
write.csv(zeros, "Results/Figure_A5/FA5-1_effects.csv", row.names = F)
ggsave("Results/Figure_A5/FA5-1_plot.pdf", width = 8, height = 8)

# get effects for counts
counts <- by_issue_mods %>%
  filter(type == "count") %>%
  ungroup() %>%
  mutate(names = issue_names,
         modelcoef = estimate,
         sign = ifelse(yhi < 0, "red", ifelse(ylo > 0, "blue", "gray"))) %>%
  dplyr::select(names, type, modelcoef, ylo, yhi, sign )

# plot
coefplot.politicized(counts) +
  ylab("Estimated Coefficient of Geopolitical Affinity")

# save
write.csv(counts, "Results/Figure_A5/FA5-2_effects.csv", row.names = F)
ggsave("Results/Figure_A5/FA5-2_plot.pdf", width = 8, height = 8)
```

# Figure 4: Effects of Issue on Probability of Recommendation Support

```{r}
# subset for relevant variables
rec_level_sub <- rec_level %>%
  dplyr::select(Response, Severity.3, ip_affinity_lagged, 
         same_region, hr_latentmean_lagged.x, hr_latentmean_lagged.y, 
         `Asylum Seekers and Refugees`:`Women's Rights`, year, To_COW, From_COW) %>%
  mutate(Severity.3 = as.factor(Severity.3),
         year = as.factor(year),
         To_COW = as.factor(To_COW),
         From_COW = as.factor(From_COW))

# model
response <- lm(Response ~ . ,
            data = rec_level_sub)
# tidy
response_tidy = tidy(response)[8:61,] %>%
  mutate(issue = issue_names) %>%
  dplyr::select(issue, estimate, std.error, p.value)

# plot
coefplot.sensitivity(response_tidy) 

# save
write.csv(response_tidy, "Results/Figure_4/F4_effects.csv", row.names = F)
ggsave("Results/Figure_4/F4_plot.pdf",  width = 8, height = 8)
```

# Figure 5: Effects of Issue Theme on Probability of Recommendation Support

```{r}
# subset for relevant variables
rec_level_sub <- rec_level %>%
  select(Response, Severity.3, ip_affinity_lagged, 
         same_region, hr_latentmean_lagged.x, hr_latentmean_lagged.y, 
         starts_with("cluster_"), year, To_COW, From_COW) %>%
  mutate(Severity.3 = as.factor(Severity.3),
         year = as.factor(year),
         To_COW = as.factor(To_COW),
         From_COW = as.factor(From_COW))

# model and tidy
response <- lm(Response ~ . ,
            data = rec_level_sub)

response_tidy = tidy(response)[8:14,] %>%
  rename("issue" = term) %>%
  mutate(issue = c("Civil-Political Rights", "Governance and Public Services", "Migration", "Physical Integrity Rights", "Racial, Ethnic, and Religious Minorities", "Socio-Economic Rights", "Protection of Vulnerable Populations"))

# plot
coefplot.sensitivity(response_tidy) + 
  theme(text = element_text(size=13))

# save
write.csv(response_tidy, "Results/Figure_5/F5_effects.csv", row.names = F)
ggsave("Results/Figure_5/F5_plot.pdf",  width = 8, height = 4)
```

# Figure 6: The Relationship between Politicization and Sensitivity

```{r}
# read coefficients from Figures 2 and 4
politicized <- read.csv("Results/Figure_2/F2_effects.csv")%>%
  dplyr::select("issue" = names, "politicization" = modelcoef) 

sensitive <- read.csv("Results/Figure_4/F4_effects.csv") %>%
  dplyr::select(issue, "sensitivity" = estimate) 

# join
issue_politicized_sensitive <- politicized %>%
  mutate(sensitivity = round(sensitive$sensitivity, 3))

# define plot
plot.rel <- function(x){
  ggplot(x, aes(x = politicization, y = -sensitivity, label = issue)) +
  geom_hline(yintercept = 0, color = "gray") +
  geom_vline(xintercept = 0, color = "gray") +
  geom_smooth(method="lm", size = .3, color = "blue") +
  geom_point() +
  xlab("") +
  ylab("") +
  geom_text_repel(aes(
    label = ifelse(abs(sensitivity) > 0.1 | abs(politicization) > 0.042, issue, '')),
    arrow = arrow(length = unit(0.02, "npc")),
    box.padding = 1,
    size = 5, color = "red") +
    theme(axis.text=element_text(size=12),
          axis.title=element_text(size=14,face="bold"))
}

# plot
plot.rel(issue_politicized_sensitive)

# save
ggsave("Results/Figure_6/F6_plot.pdf", width = 10, height = 8)

# test significance in relationship
cor.test(issue_politicized_sensitive$politicization, issue_politicized_sensitive$sensitivity)

# remove outliers - death penalty + SOGI + include race
no_out_1 <- issue_politicized_sensitive %>%
  filter(!issue %in% c("Death Penalty", "Racial Discrimination", "Sexual Orientation and Gender Identity"))

# test significance in relationship
cor.test(no_out_1$politicization, no_out_1$sensitivity)
```



# Table A1: Action and Severity Codes

```{r}
# Analysis of "Action" codes.
table_A1 <- rec_level %>%
  group_by(Action) %>%
  summarise(N = n(),
            `Percent Supported by SuR` = round(sum(Response==1)/n(), 4)*100) %>%
  mutate(`Severity Code` = c(1, 1, 3, 2, 3))

write.csv(table_A1, "Results/Table_A1.csv", row.names = F)
rm(table_A1)
```


# Table A2: Issue Codes and Frequency

```{r}
# Count recs by issue
table_A2 <- rec_level %>% 
  dplyr::select(`Asylum Seekers and Refugees`:`Women's Rights`)

table_A2 <- data.frame(Issue = names(table_A2), 
                      Counts = colSums(table_A2)) %>% 
  arrange(-Counts)

# save
write.csv(table_A2, "Results/Table_A2.csv", row.names = F)
rm(table_A2)
```

# Table A3: Number of Reviewer Dyads Issuing At Least One Recommendation about a Given Issue

```{r}
#  Table A3
table_A3 <- dyad_level %>%
  dplyr::select(`Asylum Seekers and Refugees`:`Women's Rights`) %>%
  gather(issue, n_issue, `Asylum Seekers and Refugees`:`Women's Rights`) %>%
  group_by(issue) %>%
  summarise(nonzero = sum(n_issue != 0, na.rm = T))

# save
write.csv(table_A3, "Results/Table_A3.csv", row.names = F)
rm(table_A3)
```



# Table A5: Determinants of UPR Participation

```{r}
# estimate model of participation
mod <- glm(participation ~ ip_affinity_lagged +
               hr_latentmean_lagged.x + hr_latentmean_lagged.y +
               same_region + upr_reviewed_session + HRC_member.x + 
               factor(From_COW) + factor(To_COW) + factor(year),
            family = binomial(link = "logit"),
            data = dyad_level)

# write
stargazer(mod, out = "Results/Table_A5.html",
          type = "html", 
          style = "ajps", 
          title = "Likelihood of Participating in a Given Review", 
          # se = list(logit1.se, logit2.se, logit3.se), 
          notes="Fixed SuR, Reviewer State, Year, and Issue Effects Omitted from the Table. Standard errors in parentheses.",  
          dep.var.labels = "Participation", 
          keep = 1:7,
          keep.stat = c("n"),
          covariate.labels=c("Geopolitical Affinity",
                             "Physical Integrity Rights (Reviewer)", 
                             "Physical Integrity Rights (Target)",
                             "Same Region", 
                             "UPR (Reviewer)",
                             "HRC Member (Reviewer)",
                              ""), 
          star.cutoffs = c(0.10, 0.05, .01))
```

# Table A7: Summary of Robustness Tests

```{r}
# function to read csv's of effects
read_effx <- function(path, type){
  doc <- read.csv(path)
  return(doc$sign)
}

# read F2 effects
base <- read.csv("Results/Figure_2/F2_effects.csv")%>%
  dplyr::select(names, modelcoef, "F2" = sign) 

# add robustness tests
base$FA1 <- read_effx("Results/Figure_A1/FA1_effects.csv")
base$FA2 <- read_effx("Results/Figure_A2/FA2_effects.csv")
base$FA3 <- read_effx("Results/Figure_A3/FA3_effects.csv")
base$FA4 <- read_effx("Results/Figure_A4/FA4_effects.csv")
base$FA5.1 <- read_effx("Results/Figure_A5/FA5-1_effects.csv")

# sort
base <- arrange(base, modelcoef)

# write
write.csv(base, "Results/Table_A7.csv", row.names = F)
```

# Table A8: Effect of Geopolitical Affinity on Number of Praising Recommendations (Severity Level 1)

```{r}
# estimate model
praise <- glm(n_recs ~ ip_affinity_lagged +
              upr_reviewed_session + same_region + 
              hr_latentmean_lagged.x +  hr_latentmean_lagged.y + 
              To_COW + From_COW + year,
            family = gaussian(),
            data = dyad_level_S1)

# write
stargazer(praise, out = "Results/Table_A8.html",
          type = "html", 
          style = "ajps", 
          title = "Number of Recs with Severity Level 1", 
          # se = list(logit1.se, logit2.se, logit3.se), 
          notes="Fixed SuR, Reviewer State, Year, and Issue Effects Omitted from the Table. Standard errors in parentheses.",  
          dep.var.labels = "Number of Recommendations", 
          keep = 1:6,
          keep.stat = c("n"),
          covariate.labels=c("Geopolitical Affinity", 
                             "UPR (Reviewer)", 
                             "Same Region", 
                             "Physical Integrity Rights (Reviewer)", 
                             "Physical Integrity Rights (Target)", ""), 
          star.cutoffs = c(0.10, 0.05, .01))
```

# Table A9: Reviewing Dyads Exchanging the Highest Amount of Praise 

```{r}
# which dyads exchange the most sev.1 recs?
high_praise <- dyad_level_S1 %>% 
  arrange(-n_recs) %>%
  select(From_COW, To_COW, year, n_recs) %>%
  head(14)

# save
write.csv(high_praise, "Results/Table_A9.csv", row.names = F)
```


